In this lab you will learn



library(tidyverse)
library(gridExtra) ## for plotting function grid.arrange()
library(deSolve) ## for ODE model
library(magrittr) ## for pipe symbol %<>%
library(knitr) ## for table-viewing function kable()
library(gganimate) ## for plotting animations with ggplot
theme_set(theme_bw())
theme_update(
  panel.grid = element_blank(),
  aspect.ratio = 1
)



Before we start, let’s define a concept which underlies our entire course as we study how populations grow – regulating factor.

Regulating factor

A regulating factor is any organism or environmental condition which modulates the population growth of a species. It is because of regulating factors that species do not grow exponentially in nature (Module 1). Regulating factors can be

How species interact with their regulating factors determines whether they can maintain a viable population size.


Niche theory

The basic result that two consumers cannot coexist on the same resource is one of the fundamental tenets of community ecology, and the basis for the famous competitive exclusion principle, which can be stated as

Two species occupying the same niche cannot stably coexist.

Framed this way, the principle is doing double duty as both a rule about coexistence and the definition of ecological niche. An ecological niche is whatever it is that two species cannot share and stably coexist. That sounds circular, but we break the circularity by connecting the niche to species resource needs and their impacts on those resources (Chase and Leibold 2003). As we have seen in the previous lab, two species growing on the same resource are sharing a niche. Two species sharing two resources can niche-differentiate by differing in how much they need and impact each resource. More generally, species can niche-differentiate by differing in how they interact with the different regulating factors affecting their growth. Regulating factors can be resources such as food, space, and nutrients, but can also be species enemies such as predators and parasites.


Q1. Recall from our Predator-Prey Dynamics lab the example of a green alga and a diatom cohabiting in a pond, both of which growing on the same resource, nitrate. The two species were able to coexist in the presence of a copepod, a predator which fed on both the alga and the diatom. Without the copepod, one of the species competitively excluded the other. Explain how this example relates to the competitive exclusion principle, and to the concept of niche differentiation. What are the regulating factors?

In this lab, we will explore ideas about how niche differentiation manifests in communities with many species, and how competition leads to predictable patterns in the phenotypes of species that coexist.



Figure: Species assemblages with high levels of species diversity include coral reefs, tropical rainforests, and algal communities. How do so many species coexist in the same area despite competition for the same resources, and what how does competition modulate their phenotypic biodiversity?



Fundamental vs realized niche

In an experimental study on a coffee farm in southern Mexico, ants were sampled at five different heights in the cofee plants at 0, 1, 2, 3, and 4 meters from the ground. Three species of ants were found, Pheidole protensa, Pheidole synanthropica, and Cromatogaster carinata. The proportional distribution of all the sampling sites at which a particular species was encountered is shown in the bar chart below:

Figure: Three ant species which coexist in Mexican coffee farms. Pheidole protensa by Adam Lazarus, Pheidole synanthropica by John Longino, Crematogaster carinata by C. Richart, from AntWeb



Q2. Do all species occur equally in all heights of the coffee plants? What can you hypothesize about the habitat preferences (here meaning the different areas of a coffee shrub) among the three species based on these observations?


The researchers selected two sites for removal of one of the species. In one site they removed C. carinata, and in another they removed P. synanthropica. Weeks later, they observed the following proportional distribution of the remaining two species:

Q3. Describe the difference from pre-removal observations. What can you now hypothesize about the competitive interactions among the three species?


This study suggests that competition affects species habitat distributions. Some species are competitively excluded from habitat that they are otherwise capable of colonizing, and coexist by specializing on different habitats. This is the idea behind the concepts of fundamental niche and realized niche.

We can think of this example as a case where species partition space. Niche differentiation here occurs via specialization to different environments: the lower, middle, and upper layers of vegetation on a coffee farm, each of which can be thought of as a different resource. Notice that we have three species coexisting on three resources.


Q4. Explain the difference between fundamental and realized niche, and the role played by competition in this difference.


But what about cohabiting species? Trees in a rainforest show no obvious spatial separation, yet they must niche-differentiate somehow, or they wouldn’t be stably coexisting. As we will see, competition is thought to set limits on the ability of species to occur in the same location, and to create specific patterns in the distribution of phenotypes across species that do coexist in the same habitat.


Character displacement

Consider a species with some degree of intraspecific phenotypic variation — i.e. members of the species may differ among themselves in a phenotype of interest, e.g. flower coloration. If the phenotype in question relates to competitive interactions, one might expect that the selective pressure exerted by competition will cause the phenotype of individuals living in areas where competition is potentially strong to differ from locations where competition is less of a factor. This idea led to the formulation of the concept of character displacement, according to which if two similar species have overlapping ranges, the phenotypic differences between those species will be accentuated in the areas where they do overlap, compared to areas where they occur in isolation from each other. The rationale is that individuals of species A whose phenotype is too close to those of species B will have lower fitness than individuals of species A with greater phenotypic difference from species B, since the latter would experience less competition from species B.

To recap the conditions for character displacement to occur:

  • There is phenotypic variation among individuals of the same species
  • The phenotype relates to competitive interactions
  • Competition is stronger between more similar individuals
  • The species’ range partially overlaps with that of another similar species

If all these conditions are met, we may expect individuals of species A living in the same location as species B to exhibit greater phenotpyic differences from species B than individuals of species A living in areas where species B does not occur.

Consider the yellow-and-white monkeyflower Erythranthe bicolor (formerly known as Mimulus bicolor), a bee-pollinated flower native to California. Some of its habitat overlaps with the common yellow monkeyflower, Erythrante guttata (formerly known as Mimulus guttatus), which is also pollinated by bees. As their names indicate, E. guttata is all-yellow, whereas E. bicolor has two morphs: all-yellow, and yellow-and-white.

When a bee visits a flower of one species and next visits a flower of the other species, neither flower gets pollinated. It is therefore in each plant’s interest to avoid attracting the same pollinator as the other species. One way to do that is by looking different from flowers of the other species.

This system meets all the conditions above for character displacement to occur, and indeed that is what Dena Grossenbacher and Maureen Stanton (2014) set out to investigate.


Figure: 1. The geographic distribution of Mimulus bicolor (open black circles) is restricted to the Sierra Nevada of California, whereas M. guttatus (open gray circles) occurs across western North America.From Grossenbacher and Stanton 2014; 2. Upper panel shows Mimulus guttatus and the two M. bicolor color morphs: all‐yellow and bicolor. Lower panel is a naturally occurring, mixed‐species patch with M. guttatus on the left and all‐yellow M. bicolor on the right All figures from Grossenbacher and Stanton 2014.


The researchers observed that the relative frequency of the yellow-and-white morph of E.bicolor (compared to the all-yellow morph) in areas where the other species was present or absent was as follows:



Q5: Are the results above consistent with character displacement? Explain.


In order to further connect flower color to pollination-mediated fitness, the authors set up an experiment where they planted either 11 plants of a single E. bicolor morph (monomorphic, low density treatment in the figure below), 22 plants of a single E. bicolor morph (monomorphic, high density), or 11 plants of one E. bicolor morph intermixed with 11 E. guttata plants (heterospecific). In each treatment, they measured the number of seeds produced per each flower (seeds develop after bees pollinate the flowers). They report the following results:


Figure: Seed set of E. bicolor yellow-and-white morphs (open bars) and all‐yellow morphs (gray bars) under three conditions with and without E. guttata: arrays containing either 11 plants of one E. bicolor morph, 22 plants of one E. bicolor morph, or 11 plants of one M. bicolor morph intermixed with 11 E. guttata plants.


Q6: Focusing on the difference between the third treatment (where E. guttata is present) from the other two treatments, explain how those results connect to the idea that flower color affects pollinator-mediated fitness in E. bicolor.


Classic limiting similarity

The notion that coexisting species must occupy different niches has an intuitive implication: coexisting species will differ from each other in some non-random way, and as such will be more different than expected by chance. This is the logical basis for the principle of limiting similarity, an influential idea in community ecology that has guided research on coexistence for decades.

Here we will explore the origins of the limiting similarity principle by numerically simulating the Lotka-Volterra (LV) model in a community with three species. We will be following in the footsteps of Robert MacArthur and Richard Levins (1969), who demonstrated that the Lotka-Volterra model predicts a limit for the similarity among three coexisting species. Specifically, they showed that if a community has two resident species, then a third one with an intermediate phenotype can only invade and coexist with the other two if its niche overlap with the residents was sufficiently low, a condition which required that the residents had sufficiently different phenotypes (hence the phrase limiting similarity).

Let’s remind ourselves of the Lotka-Volterra competition equations generalized for an arbitrary number of species:

\[ \frac{dN_i}{dt} = r_iN_i\left(1-\frac{1}{K_i}\sum_j \alpha_{ij}N_j\right) \tag{1} \] which you can compare with the LV model for two species (Equation 3 in the previous lab). The parameters are the intrinsic growth rate \(r_i\) and carrying capacity \(K_i\) of each species, and the competition coefficients \(\alpha_{ij}\) which specify how much species \(j\) regulates species \(i\) relative to how much species \(i\) regulates itself (by definition, \(\alpha_{ii} = 1\)).

The R code below provides functions for implementing the model and plotting its outcome:

LV_Model = 
  function(
    initial_N,
    intrinsic_growth_rate,
    carrying_capacity,
    community_matrix,
    final_time,
    time_step
  ){
      LV = 
        function(t, state, parameters){
           with(as.list(parameters), {
             N = state
             dNdt = as.numeric(r * N * (1 - A %*% N / K))
             list(dNdt)
           })
          }
      
      times = seq(0, final_time, by = time_step)
      
      parameters = 
        list(
          r = intrinsic_growth_rate,
          K = carrying_capacity,
          A = community_matrix
        )
      
      state = initial_N
      
      out = ode(y = state, times = times, func = LV, parms = parameters)
      
      return(
        list(
          parameters = 
            list(
              intrinsic_growth_rate, 
              carrying_capacity, 
              community_matrix
            ),
          initial_conditions = initial_N,
          state = out
        )
      )
   }

# Plotting the model outcome
Plot_LV_timeseries = 
  function(model){
    as.data.frame(model$state) |>
    pivot_longer(-time, names_to = 'species', values_to = 'abundance') |>
    ggplot(aes(time, abundance, group = species, color = species)) +
    geom_line(size = 1) +
    expand_limits(y = 0)
  }

Plot_LV_stemplot = 
  function(trait, model){
    tibble(trait = trait, abundance = round(model$state[nrow(model$state), -1])) |>
    ggplot() +
    geom_segment(aes(trait, abundance, xend = trait, yend = abundance - abundance)) +
    geom_point(aes(trait, abundance))
  }


At this stage we must recognize that while the LV model tells us that species \(j\) affects species \(i\) via the competition coefficient \(\alpha_{ij}\), it has nothing to say about what those coefficients should be. We now turn to this task.

The basic idea is to write down a model that specifies how much species compete based on their respective phenotypes. We will do so by assuming a link between a species phenotype (ie its trait) and its interaction with the different resources. For the sake of example, let’s think back on Darwin’s finches, famous for the variety of their beak shape/size and the corresponding specialization to different diets. Let’s simplify things by assuming that all finches eat seeds of different sizes that they find on the islands, but each species has its own preferences for a certain seed size based on the shape/size of its beak.

Suppose we line up our finch species on a trait axis based on their beak size, from the species with the smallest beak on one extreme of the axis to the species with the largest beak on the other extreme (ie the trait here is beak size). Making the very reasonable assumption that birds with similar beak sizes have similar seed size preferences, we presume that similarity in beak size (trait) implies similarity in diet (resource use), and therefore stronger niche overlap, and therefore more intense competition. In other words, if we were to plot the competition coefficient between a pair of species against the difference in beak size, the curve would peak at zero trait difference and taper off at higher trait difference.

Figure: Top. Four of Darwin’s finches, showcasing the gradation of beak shape/size. Illustration by John Gould (1804-1881). Left. The beak shape corresponds to the species diet: small-beaked birds preferentially eat small seeds, while large-beaked birds prefer large seeds. The overlap in diet leads to competition. Right. Because similarity in beak shape reflects similarity in diet, birds with more similar beak shape will compete more intensely than birds with very different beak shapes.


Q7. Following the logic of the paragraph above and based on the shape of the competition curve in the bottom-right plot above, which finch species is the biggest competitor (“worst enemy”) of Geospiza fortis? (consider all species, including G. fortis)


Q8. Recalling what you know about how the prospects for two-species coexistence relate to intraspecific vis-à-vis interspecific competition, would you say the scenario above is one where coexistence is possible in principle? Explain your answer.


One way to implement the ideas above is to set the competition coefficients as follows:

\[ \alpha_{ij} = e^{-\left(\frac{x_i - x_j}{w}\right)^4} \tag{2} \] In words, this expression is saying that a species \(j\) with trait (beak size) \(x_j\) competitively affects species \(i\) with trait \(x_i\) via a function that decreases exponentially with the trait difference \(x_i - x_j\), scaled by the parameter \(w\). The larger the \(w\), the more the species compete. MacArthur showed that an expression similar to this one arises as a measure of the overlap in the resource utilization curves depicted in the figure above for our finches-and-seeds example.

Figure: Visual representation of the competition coefficient \(\alpha_{ij}\) as a function of the trait difference between species \(i\) and \(j\) defined in Equation (2). Species with similar traits compete more strongly than species with dissimilar traits.


Q9. According to the expression above, what is the value of \(\alpha_{ii}\) (i.e. the intraspecific competition coefficient)?


We are now ready to implement our model of competition based on traits. In order to focus on the effect of our assumption that competition depends on trait similarity, let’s assume all species have the same intrinsic growth rates and carrying capacities regardless of their beak size.


Q10. If a finch species is alone on the island, its population will grow to its carrying capacity \(K\). Assuming \(K\) is related to the amount of resource on the island, what does the assumption of equal \(K\)’s across finches correspond to in terms of the availability of seeds of different sizes on the island?

Let’s see what happens when we start a community with two resident species at their carrying capacity and introduce a third species with 10 individuals. We assume the residents have trait values -0.5 and 0.5, and the invader has intermediate trait value 0. We run the model for 500 time units (e.g. years). Let’s set the scaling parameter \(w = 0.2\).

number_of_species = 3
scale = .2
trait = c(-.5, 0, .5)
K = 1000
N0 = c(K, 10, K)
r = .1
final_year = 500

trait_differences = as.matrix(dist(trait))
community_matrix = exp(-(trait_differences / scale) ^ 4)

model = 
  LV_Model(
    initial_N = N0, 
    intrinsic_growth_rate = c(r, r, r), 
    carrying_capacity = c(K, K, K), 
    community_matrix = community_matrix, 
    final_time = final_year, 
    time_step = 1
  )

plot_timeseries = 
  Plot_LV_timeseries(model) + 
  theme(legend.position = 'none') +
  theme(aspect.ratio = 1)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot_stemplot = 
  Plot_LV_stemplot(trait = trait, model = model) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  theme(aspect.ratio = 1) +
  xlim(c(-.5, .5))

gridExtra::grid.arrange(plot_timeseries, plot_stemplot, nrow = 1)

The first plot shows the abundances of each of the three species over time (blue are the residents, green is the invader), and the second plot shows the final abundance of each species against their trait values.

As you can see, the third species has no issues invading and coexisting with the residents. In fact, the residents don’t even seem to respond to the invader, as their populations are not affected by the invasion.

Let’s now assume the residents are more similar than before, placing them at trait values -0.2 and 0.2, and keep everything else the same:

number_of_species = 3
scale = .2
trait = c(-.2, 0, .2)
K = 1000
N0 = c(K, 10, K)
r = .1
final_year = 500

trait_differences = as.matrix(dist(trait))
community_matrix = exp(-(trait_differences / scale) ^ 4)

model = 
  LV_Model(
    initial_N = N0, 
    intrinsic_growth_rate = c(r, r, r), 
    carrying_capacity = c(K, K, K), 
    community_matrix = community_matrix, 
    final_time = final_year, 
    time_step = 1
  )

plot_timeseries = 
  Plot_LV_timeseries(model) + 
  theme(legend.position = 'none') +
  theme(aspect.ratio = 1)

plot_stemplot = 
  Plot_LV_stemplot(trait = trait, model = model) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  ) +
  theme(aspect.ratio = 1) +
  xlim(c(-.5, .5))

gridExtra::grid.arrange(plot_timeseries, plot_stemplot, nrow = 1)

Notice that now the invader is still able to invade, but its population never reaches its carrying capacity (ie the population level it would achieve if growing alone). The resident populations now decrease in response to the invasion. But ultimately all three species still coexist.


Q11. Try the experiment again, now with the residents at -0.15 and 0.15. What happens? Show your plots and describe the outcome.


We conclude that the trait difference between the invader and the residents strongly affects the outcome. Namely, if the two residents are too close on the trait axis (i.e. beyond a certain limit to similarity), coexistence with a third intermediate species is not possible. In fact, MacArthur and Levins showed that if that difference \(d\) is much smaller than the scaling parameter \(w\), ie if \(d \ll w\), the invader cannot coexist.

Note: Notice that in this trait-based Lotka-Volterra model, the number of resources (types of seeds) vastly exceeds the number of consumers (finch species), and could be, in fact, infinite (i.e. there could be continuous variation in seed size). And still, we see that coexistence is not always possible. This reflects the fact that while having at least as many resources as consumers is necessary for coexistence, it is not sufficient.

What is the biological meaning of the scaling parameter \(w\)? MacArthur showed that Equation (2) above for the competition coefficients arises as a measure of the overlap in the resource utilization curves of the consumers (bottom left panel in the figure above). The scaling parameter relates to the width of those curves. The wider those curves, the wider the birds’ diets, and the higher the scaling parameter \(w\). In turn, a lower \(w\) means that the residents are more specialized on their preferred seed size. More specialized birds can be closer together on the trait axis—ie have more similar beak sizes—and still coexist with the invader. In other words, the scaling parameter reflects the degree of resource specialization, and sets the relationship between trait difference and niche overlap. The more specialized (smaller \(w\)), the less niche overlap for the same trait difference.

This \(d_\textrm{min} \approx w\) result was the first quantitative formulation of the limiting similarity principle, and provided a testable theoretical prediction for coexistence in nature. (Note: we simplified things by assuming all species have the same parameters, but one can prove limits to similarity for the more general case of different carrying capacities and intrinsic growth rates. Feel free to run the model using different \(K\)’s for the three species.)

However, empirical evidence for this rule was limited, and eventually the \(d_\textrm{min} \approx w\) rule was relaxed to a less restrictive prediction that coexisting species should be more different from each other than a null hypothesis of random trait differences. But even that prediction was met with mixed evidence. In fact, for every study that found greater-than-chance differences between species, one or more studies found the opposite! (For more details see D’Andrea and Ostling 2016.)

In the next section we will see that one of the reasons for this mixed evidence is that the predictions for the traits of coexisting species change when we consider high-diversity systems subject to immigration from outside and demographic stochasticity.


Generalized limiting similarity: Trait clustering

It is commonly said that one square mile of rainforest may contain as many plant species as all of Europe. Despite all this diversity, all these plants are competing for the same handful of resources: light, water, and nutrients. This seems to fly in the face of our coexistence rules relating the number of species to resources. How do we explain this coexistence? And what role does competition play in structuring this prodigious biodiversity? Here we will use our model to explore coexistence in high-diversity communities, with tropical rainforests in mind.

What happens when we run the model with more than three species?

Let’s try 100 species. For simplicity we will keep carrying capacities and intrinsic growth rates the same for all species.

number_of_species = 100
scale = .2
trait = seq(-.5, .5, length = number_of_species)
K = 1000
N0 = rep(K, number_of_species)
r = .1
final_year = 5000

trait_differences = as.matrix(dist(trait))
community_matrix = exp(-(trait_differences / scale) ^ 4)

model = 
  LV_Model(
    initial_N = N0, 
    intrinsic_growth_rate = rep(r, number_of_species), 
    carrying_capacity = rep(K, number_of_species), 
    community_matrix = community_matrix, 
    final_time = final_year, 
    time_step = 1
  )


Since a plot with 100 abundance-by-time curves is not particularly illuminating, we will instead plot species by their trait values and animate their abundances over time.

We notice several interesting things from our result above.

  1. Even though we start out with 100 species, many species quickly go extinct.

  2. The remaining species quickly sort into groups.

  3. Those groups are well separated from each other. In fact, each group is maximally separated from its neighbors on either side.

  4. Eventually only one species per group will survive.

  5. However, once only a few species remain, this species sorting process slows down considerably.


Q12. Recall that in this model, trait similarity between species determines the intensity of competition. With this in mind, how do you explain the fact that the groups are maximally separated?


We can think of this pattern as the extension of the idea of limiting similarity to a multispecies community. The focus here is not on a specific limit to the similarity of coexisting species, but on the idea that trait separation between coexisting species will tend to be higher than expected by chance.

The result above also introduces another way to think about ecological niches. When the community finally reaches its equilibrium where species abundances are no longer changing, there will be five species left coexisting in stable equilibrium. Since ecological niches are conceptually associated with stable coexistence, one could say this community has five niches, and therefore five species can stably coexist.

However, another prominent feature of the result above is the timing of competitive exclusion. While the niches quickly sort themselves out with exclusion zones between them, competitive exclusion inside each niche is slow. If we take a snapshot of the community before it reaches its final stage, what we see is well separated niches, each being shared by mutually similar species.

In sum, we have learned that

  1. In equilibrium, there can only be one species per niche. One species, one niche.

  2. However, before this equilibrium is finally reached, each niche is shared by mutually similar species among which competitive exclusion is slow.


Q13. Set the scaling parameter to 0.1 and rerun the model (use the code below and substitute the ??). How many niches do we have now?

number_of_species = 100
scale = ??
trait = seq(-.5, .5, length = number_of_species)
K = 1000
N0 = rep(K, number_of_species)
r = .1
final_year = 5000

trait_differences = as.matrix(dist(trait))
community_matrix = exp(-(trait_differences / scale) ^ 4)

model = 
  LV_Model(
    initial_N = N0, 
    intrinsic_growth_rate = rep(r, number_of_species), 
    carrying_capacity = rep(K, number_of_species), 
    community_matrix = community_matrix, 
    final_time = final_year, 
    time_step = 1
  )

Plot_LV_stemplot(trait = trait, model = model)


As previously mentioned, the scaling parameter \(w\) reflects how quickly interspecific competition drops as trait separation increases, which depends on the degree of specialization of the consumers. Here we can see that \(w\) is related to the width of the niches measured as a range of trait values. The smaller the \(w\), the narrower the niches, and therefore the more niches can be packed within the same trait range.


Trait clustering under stochastic dynamics

As said many times in this course, the Lotka-Volterra competition model gives us valuable insight into the rules of coexistence, but we must keep in mind that it is a simplistic model. Before we theorize that its lessons are relevant for the real world, we must make sure that its predictions are robust to real-life complexity left out of the model.

In particular, recent developments in community ecology have emphasized the importance of stochastic (ie random) forces as drivers of ecological dynamics. Two of the most important stochastic forces are dispersal and ecological drift. In simple terms, we could say that real life communities are not isolated, and are subject to random events. As such, we must account for immigration, and we must think of birth and death as probabilistic rather than deterministic processes.

The R code below rewrites our Lotka-Volterra model under this stochastic framework.

Stochastic_LV_Model = 
  function(
    trait,
    initial_N,
    intrinsic_growth_rate,
    carrying_capacity,
    immigration_pressure,
    community_matrix,
    final_time,
    turnover_rate
  ){
    
    S = length(trait)
    r = intrinsic_growth_rate
    K = carrying_capacity
    A = community_matrix
    N = initial_N
    m = immigration_pressure
    
    N_record = N
    
    time = 0
    
    while(time < final_time){
      
      time = time + 1
      deaths = 
        table(
          sample(
            S, 
            replace = TRUE, 
            size = turnover_rate * sum(N), 
            prob = r * N
          )  
        )
        
      death_ids = as.numeric(names(deaths))
      
      births = 
        table(
          sample(
            S, 
            replace = TRUE, 
            size = (1 - m) * turnover_rate * sum(N), 
            prob = pmax(0, r * N * (2 - 1/K * as.numeric(A %*% N)))
          )
        )
        
      birth_ids = as.numeric(names(births))
      
      immigrants = 
        table(
          sample(
            S,
            replace = TRUE,
            size = m * turnover_rate * sum(N),
          )
        )
      
      immigrant_ids = as.numeric(names(immigrants))
      
      N[death_ids] = pmax(0, N[death_ids] - deaths)
      N[birth_ids] = N[birth_ids] + births
      N[immigrant_ids] = N[immigrant_ids] + immigrants
      
      N_record = rbind(N_record, N)
    }
    
    return(
      list(
        parms = list(r, K, A, m, N0 = initial_N, final_time),
        state = cbind(0:final_time, N_record)
      )
    )
    
}


The model works as follows: each year, a percentage of the trees die and are replaced by recruitment of new seedlings. (Since all deaths are replaced, the total number of trees is fixed. Relaxing this simplifying assumption does not change results.) Each tree has the same probability of death as all other trees. However, the per capita probability of recruitment is affected by competition based on trait similarity between species, just like in the deterministic LV model above. In addition, a proportion of the recruits are seedlings dispersing into the community from outside—immigrants. All species have equal probabilities of immigrating every year.

The goal here is to see whether the species groups, or clusters, that emerged from the LV dynamics above still appear under this more realistic scenario.

To keep with the spirit of realism, we will use parameter values consistent with a real-world species assemblage. Specifically, we will focus on Barro Colorado Island, Panama, which we encountered in Module 6 while discussing spatial patterns. As you may recall, a 1,000m by 500m plot (50 hectares) at the center of the island has been censused at regular intervals for decades. The plot has about 20,000 adult trees spanning 200 species for which both trait and abundance data are available. We therefore start our community with 20,000 trees from 200 species (100 trees each), and give them all a carrying capacity \(K=20,000\), meaning that if the species were alone it would reach a population of 20,000 trees. Also, it has been estimated that 1% of recruiting trees on the 50-hectare plot are immigrants. Therefore we set the immigration pressure parameter \(m = 0.01\).

The scaling parameter \(w\) in the competition coefficients \(\alpha_{ij}\) is less easy to ground on observations, since those are not easily measurable. As above, we will standardize the trait values to ranging between -0.5 and 0.5, and we will arbitrarily set the scaling parameter \(w = 0.2\).

Figure: Barro Colorado Island, by Chris Ziegler



number_of_species = 200
scale = .2
trait = seq(-.5, .5, length = number_of_species)
K = 20000
N0 = 100
r = .1
m = .01
final_year = 5000

trait_differences = as.matrix(dist(trait))
trait_differences = pmin(trait_differences, 1 - trait_differences)
community_matrix = exp(-(trait_differences / scale) ^ 4)


We are now ready to run the model:

set.seed(0)

model_stochastic = 
  Stochastic_LV_Model(
    trait = trait,
    initial_N = rep(N0, number_of_species), 
    intrinsic_growth_rate = rep(r, number_of_species), 
    carrying_capacity = rep(K, number_of_species), 
    immigration_pressure = m,
    community_matrix = community_matrix, 
    final_time = final_year, 
    turnover_rate = 0.1
  )


dtf_stochastic = 
  model_stochastic$state |> 
  as_tibble() |> 
  pivot_longer(-V1, names_to = 'species', values_to = 'abundance') |> 
  rename(year = V1) |> 
  mutate(trait = trait[match(species, unique(species))])


Like before, we can plot the species by their trait value and animate their abundances:


We notice similarities and differences from the fully deterministic LV model.

  1. The clusters are still visible.

  2. The clusters are noisy.

  3. The clusters are no longer transient, but rather a permanent feature in the community.


Q14. Why are the clusters noisy rather than neat and predictable as before (no need to plot anything, just explain the noise in the animation above)?


Q15. In this simulation, we also incorporated immigration of species. Consider for a minute how immigration should affect the dominance or extinction of species within clusters. Use that thought to explain why the clusters are now permanent rather than transient (i.e. temporary) as before.

If we can census the community a large number of times, we can calculate and plot the time-averaged abundance of each species:

plot_stem_means = 
  dtf_stochastic |> 
  filter(year > 2500) |>
  group_by(trait, species) |>
  summarize(abundance = mean(abundance), .groups = 'drop') |>
  ggplot() +
  geom_segment(aes(x = trait, y = abundance, xend = trait, yend = abundance - abundance)) +
  geom_point(aes(trait, abundance))

plot_stem_means


Notice how all three dynamic forces—competition, immigration, drift—leave their signature on the community: Competition creates the niches to which species are sorted, which appear as species clusters; immigration maintains the clusters permanently; drift adds a random element to the cluster shapes.

This finding was the theoretical underpinning of D’Andrea et al. 2020, who found evidence of this type of clustering on Barro Colorado Island.

Figure: The 50-hectare Forest Dynamics Plot in Barro Colorado Island displays species clustering by maximum height. Alternating colors identify the clusters.



The authors found that BCI trees are segregated into four groups by maximum height (ie the maximum height that a tree of a certain species can achieve): shrubs, understory treelets, midstory trees, and canopy trees. This result was interpreted as evidence that competition for light is an important force shaping the vertical structure of the forest. They concluded that BCI has four “light niches”.

Figure: Visual representation of the forest on Barro Colorado Island, highlighting the four height-based groups: shrubs (61 species; 8,071 adult trees), understory treelets (71; 7,308) midstory trees (61; 2,957), and canopy trees (64; 2,836). From D’Andrea et al. 2020



Q16: According to the sunfleck model (Terborgh 1985), light incidence below the canopy in forests does not decrease steadily towards the ground, but rather has local peaks at intermediate heights, corresponding to the convergence points of sun rays coming in at an angle. In the case of tropical forests, there are three such sub-canopy light peaks. How might the results of D’Andrea et al 2020 relate to the sunfleck model?


Figure: Schematic model of light penetration into an eastern oak-dogwood forest. Diagonal lines represent sun rays shining on the forest at specific angles at different times of day. Heights below the intersection closest to the ground (~7m from the ground) receive the same total amount of sunlight. However, daily light incidence is gappy above that spot, peaking at the intersections shown. This creates the possibility of light-niches at specific heights above the ground. Due to differences in incidence angles and canopy shapes, the model predicts 3 sub-canopy intersections in tropical forests, rather than the 2 shown here for a temperate forest. From Terborgh 1985.


Q17. According to the sunfleck model, temperate forests have not three but two subcanopy light peaks, due to differences in canopy shape and sunlight angle of incidence compared to tropical forests. How many height clusters would you expect to see by applying the methodology of D’Andrea et al 2020 to Harvard Forest in Massachusetts?


Figure: Harvard Forest in Petersham, Massachusetts, is an example of a temperate forest. Photo by Maggie Janik.


Q18. Harvard Forest contains about 25 tree species. By comparison, the BCI plot has 257. What challenges might you expect a researcher would face upon analyzing Harvard Forest for species clusters? Hint: Cluster analysis requires a large number of members per cluster, otherwise it becomes difficult to reject the null model.

—